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Abstract 

The  unbiased  FIR  filter  is  investigated  to  estimate  the  time  interval  error  (TIE)  K-D  polynomial 
model  of  a  local  clock  in  GPS-based  timekeeping  in  the  presence  of  uniformly  distributed  sawtooth  noise. 
An  estimation  algorithm  is  proposed  and  applied  to  the  GPS-based  TIE  measurements  of  a  crystal  clock 
without  using  the  sawtooth  correction.  Based  upon  this,  we  show  that  the  TIE  estimates  fit  actual  values 
with  an  uncertainty  of  the  GPS  time.  It  is  also  demonstrated  that  estimates  of  the  fractional  frequency 
offset  fit  the  measurements  with  the  frequency  shifts  present  in  the  reference  rubidium  source  and  the  1 
PPS  signal  of  a  GPS  receiver  used. 


INTRODUCTION 


Fast  and  accurate  estimation  and  steering  of  a  local  clock  performance  is  still  a  key  problem  in  GPS-based 
timekeeping.  Here,  the  time  interval  error  (TIE)  between  the  GPS  time  and  the  local  clock  time  is  measured 
in  the  presence  of  noise  induced  by  the  GPS  timing  receiver.  The  TIE  is  then  estimated  and  the  feedback 
system  causes  the  local  clock  to  track  the  GPS  time.  The  standard  deviation  of  the  noise  using  commercially 
available  receivers  is  about  30  ns,  can  reach  10-20  ns  [1] ,  and  may  be  improved  by  removal  of  systematic  errors 
to  no  less  than  3-5  ns  [1,2].  Having  such  a  large  noise,  the  measured  data  cannot  be  used  straightforwardly 
for  locking,  and  a  TIE  tracking  filter  is  applied.  To  obtain  filtering  in  an  optimum  way,  a  TIE  model  of  a 
local  clock  must  be  known  for  the  filter  memory.  Such  a  model  was  proposed  in  [4]  as  a  third  degree  (3-D) 
Taylor  polynomial  and  is  now  basic  in  timekeeping,  being  practically  proved.  In  the  discrete  time,  the  model 
may  be  written  as 


x(n) 


xo  +  yo^n  +  yA2n2  +  w(n,r), 


(1) 


where  n  =  0, 1, ...;  A  =  tn  —  f„_i  is  a  sample  time;  tn  is  a  discrete  time;  Xq  is  an  initial  time  error;  yo  is  an 
initial  fractional  frequency  offset  of  a  local  clock  from  the  reference  frequency;  D  is  an  initial  linear  fractional 
frequency  drift  rate;  and  w(n,  r)  is  a  random  component  caused  by  the  oscillator  colored  Gaussian  noise  and 
environmental  influences.  In  GPS-based  measurements,  (1)  is  observed  via  the  mixture 


£(n)  =  x(n)  +  v(n),  (2) 

in  which  v(n)  is  a  noisy  component  induced  at  receiver  (the  noise  of  a  measurement  set  is  usually  negligible). 
It  is  also  typical  for  GPS-based  measurements  that  the  noise  mean  square  values  relate  as  ( w2(n,T ))  << 
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( v2(n )).  Therefore,  w(n,r)  may  be  neglected  in  the  averaging  procedure  (it  cannot  be  discarded  in  the 
Kalman  filter). 

In  modern  receivers,  such  as  the  Motorola  M12+  (see  [3]),  and  SynPaQ  III  GPS  Sensor,  a  random  variable 
v[n)  is  uniformly  distributed  owing  to  a  principle  of  the  1  PPS  (one  pulse  per  second)  signal  formation. 
Figure  la  shows  a  typical  sawtooth  structure  of  v(n)  in  a  short  time;  that  is  the  modulo  2 amax  Brownian 
TIE  associated  with  a  phase  of  a  receiver  local  oscillator,  where  amax  is  the  error  bound  (in  SynPaQ  III,  it 
is  about  50  ns).  In  long-term  measurements,  v(n)  exhibits  nonstationary  excursions  (Figure  16)  caused  by 
uncertainties  of  a  GPS  time  at  receiver  (see  [5])  with  different  satellite  sets  in  a  view.  If  a  single-channel 
receiver  is  used,  the  excursions  have  a  day’s  periodicity  [5] .  With  multichannel  receivers,  they  are  typically 
nonperiodic  (Figure  16).  In  the  latter  case,  a  long-term  noise  histogram  evolves,  by  increasing  the  data,  from 
uniform  to  normal  and  vn  may  be  approximated  by  a  mean-zero,  ( v(n ))  =  0,  stationary  Gaussian  noise.  If 
the  sawtooth  correction  is  used,  then  the  noise  v(n)  becomes  near  Gaussian  in  a  short  time. 

To  estimate  optimally  Xq,  yo,  and  D  for  different  clocks  (atomic  and  crystals),  assuming  the  Gaussian 
nature  of  v(n)  with  a  known  variance  <r]j,  several  filtering  algorithms  have  been  examined  over  a  couple  of 
decades.  For  the  state  space  equivalent  of  (1),  the  problem  was  formulated  by  Allan  and  Barnes  in  [6]  to 
apply  Kalman  filtering.  Several  solutions  were  then  given  in  [7-9].  Thereafter,  various  Kalman  algorithms 
were  examined  by  the  authors  in  [10-14]  for  different  estimation  purposes.  These  and  other  applications 
of  Kalman  filters  to  time  scales  were  recently  outlined  in  [15].  The  problem  with  the  Kalman  filter  arises 
when  noise  is  non-Gaussian,  as  in  Figure  1,  for  example.  To  overcome  this,  advanced  Kalman  algorithms 
were  proposed  in  [16,17],  being,  however,  not  yet  adopted  for  time  scales.  Alternatively,  finite  impulse 
response  (FIR)  filters  are  used  that  allow  for  noise  of  an  arbitrary  distribution.  In  contrast  to  the  infinite 
impulse  response  (HR)  structures  (including  Kalman  filters),  FIR  structures  have  inherent  properties,  such 
as  a  bounded  input/bounded  output  (BIBO)  stability  and  robustness  against  temporary  model  uncertainties 
and  round-off  errors  [20].  They  may  be  used  independently  or  combined  with  Kalman  filters  [12,19].  A 
general  optimal  FIR  filtering  algorithm  with  embedded  unbiasedness  for  state  space  models  was  recently 
proposed  in  [21].  Especially  for  GPS-based  timekeeping  and  a  linear  TIE  model,  an  unbiased  FIR  filter  was 
designed  and  studied  in  [22].  In  this  paper,  we  investigate  a  new  unbiased  FIR  filter  and  filtering  algorithm 
intended  for  real-time  estimation  of  the  TIE  AT-D  polynomial  model  in  the  presence  of  noise  of  an  arbitrary 
symmetric  distribution. 


THE  TIE  POLYNOMIAL  MODEL 


In  view  of  (ut^)  «  (r^),  the  noise-free  and  time-invariant  TIE  model  projects  ahead  on  a  horizon  of  N 
points  from  the  start  point  n  =  0  with  the  K-D  Taylor  polynomial 


I< 

xi(n)  = 

p—0 


A  Pnp 

V ! 


cci(O)  +  £2(0)A n  - 


x3(n) 


A  2n2  + 


X4  (n) 


A6n6.. 


(3) 


By  extending  the  time  derivatives  of  x\(n)  to  the  Taylor  series,  the  clock  model  and  its  observation  (2) 
become,  respectively, 

A(n)  =  B(n)A(0),  (4) 

((n)  =  CA(n)+v(n),  (5) 

where  A (n)  =  [x\{n)x2{n) ...XK+i{n)]T  is  a  vector  [( K  +  l)xl]  of  clock  functions  with  the  components 
approximately  calculated  for  l  >  1,  l  €  [0,  A'],  by 

xi(n)  =  ^[xi-ifo)  -  xi-i(n  -  1)],  (6) 
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and  a  transition  matrix  \{K  +  1)  x  (K  +  1)]  is 


'  1 

An 

A2n2/2  . 

.  (An)K/K\ 

0 

1 

An 

.  {An)K~l /{K  - 

-1)! 

B(n)  = 

0 

0 

1 

.  (A n)K~2/{I<  - 

-2)! 

_  0 

0 

0 

.  1 

The  most  common  situation  that  may  be  assumed  in  timekeeping  is  when  all  or  several  clock  states  are 
observable,  by  (5),  via  M  (independent  or  dependent)  measurements  in  presence  of  correlated  or  uncorrelated 
noises.  Then  the  observation  vector  is  £(n)  =  [£i(u)£2(u)...<!;m(u)]T  and 


cii  0  ...  0  ...  0 

0  C22  0  ...  0 

.  0 

0  0  ...  cmm  •  0 


(8) 


is  the  measurement  matrix  [M  x  ( AT  +  1)]  with,  typically,  unit  components.  Finally,  in  view  of  Figure  1, 
the  noise  vector  v(n)  =  [ui(?i)u2(ti)...Um(u)]t  contains  correlated  or  uncorrelated  components  that  are  not 
obligatory  Gaussian.  In  the  following,  we  will  assume  that  the  TIE  model  (4)  has  K  +  1  states  and  is 
observable  via  (5)  in  presence  of  sawtooth  noise  (Figure  la)  of  unknown  structure  (Figure  16)  caused  by 
GPS  time  uncertainty.  We  will  also  assume  that  M  =  K  +  1,  C  is  an  identity  matrix,  and  the  noise  v(n)  is 
mean  zero,  (v(n))  =  0,  and  symmetrically  distributed  with  known  covariance  (v 1  (n)v(n)). 


AN  UNBIASED  FIR  FILTER 


Utilizing  N  points  of  the  nearest  past,  the  FIR  estimate  A (n)  =  [xi(n)x2(n)...Xx+i(n)]T  at  n-th  point  is 
given  by,  for  M  =  K  +  1, 

JV-l 

A(n)=  £H(*)£(n-i)  (9) 

i= 0 

=  q(n)r,  (10) 

=  d(n)r  +  r(n)r,  (11) 

where  the  matrix  [(A'  +  1)  x  (AT  +  1]  of  unknown  FIRs  is 


H(*)  = 


hK(i)  0 
0  hK-i(i ) 


0 


0 


in  which  the  Z-tli  FIR  has  inherent  properties:  hi(i)  = 


0 
0 

...  h0(i) 

hi  (i) ,  0  <  *  <  N  —  1 
0,  otherwise 


(12) 


and  =  1-  The 


estimates  vector  and  its  deterministic  and  random  constituents,  all  of  [(AT  +  1)  x  N],  are,  respectively, 


q(n) 
d(n)  = 
r(n) 


=  [H(0)£(n) 

H(l)^(n  —  1)  .. 

.  H(A-l)$(n-  JV+1)] 

(13) 

[H(0)CA(n) 

H(l)CA(n  —  1)  . 

...  H(N  —  l)CA(n  —  N  +  1)], 

(14) 

=  [H(0  )v(n) 

H(l)v(n  —  1)  .. 

.  H(7V  —  l)v(n  —  N  +  1)], 

(15) 
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and  an  auxiliary  unit  matrix  (N  x  1)  is  T  =  [1  1  ...1]T.  The  mean  square  error  (MSE)  of  A (n)  is 

J(n)  =  ([A (n)  -  A(n)f[A(n)  -  A(n)]>  =  [A(»)  -  d(n)T]T[A(n)  -  d(n)T]  +  <[r(n)r]T[r(n)r]),  (16) 

in  which  the  estimate  bias  and  variance  are,  respectively, 

bias[A(n)]  =  A  (n)  —  d(n)r,  (17) 

var[A(n)]  =  ([r(n)r]T[r(n)r]).  (18) 


Generic  Coefficients  for  the  FIR  of  an  Unbiased  Filter 

The  necessary  and  sufficient  condition  for  the  unbiased  estimate  follows  straightforwardly  from  (17);  that  is, 


A(n)  =  d(n)T, 


providing  the  rule  to  derive  FIRs  for  the  clock  states,  namely: 


xi{n) 

'  W^-Afin) 

X2  (n) 

= 

W^_1A2(n) 

_  xK+i  (n) 

_  W^A_fc+i(n)  _ 

where 


W,  =  [/q(0)/q(l).../q(AT-l)]T, 


Aic+i-;(n) 


xK+i -i(n) 
xK+i -i(n  -  1) 

xK+i-i(n  —  N  +  1) 


For  the  clock  (K  +  1  —  l)- th  state,  (20)  thus  yields  a  relation 


xK+i -i(n)  =WfXK+1-i(n). 


(19) 

(20) 

(21) 

(22) 

(23) 


By  invoking  (6)  and  then  expressing  the  components  in  (20)  with  the  Z-D  polynomial,  by  (3),  we  arrive  at 
the  unbiasedness  (or  deadbeat)  constraint 

FW;  =  G,  (24) 

in  which  G  =  [1  0  ...  0]T  is  of  [(l  +  1)  x  1]  and 


'  1 

i 

i 

...  1 

0 

i 

2 

...  IV -1 

0 

i 

22 

...  ( AT-1 )2 

o  : 

i 

2l 

...  {N-l)1 

We  notice  that  the  constraint  (24)  is  inherent  for  any  other  unbiased  estimators,  e.g.,  for  the  minimum 
variance  unbiased  (MVU)  and  best  linear  unbiased  estimator  (BLUE).  Now,  the  components  in  (21)  may 
also  be  substituted  by  the  Z-D  polynomial  (that  is  what  Kalman  claims  for  optimal  filtering) 


i 

hi(i)  =  T^ajiiP, 

3=0 


(26) 
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where  aji  are  still  unknown  coefficients.  Embedded  (26),  the  constraint  (24)  becomes 


DA  =  G, 


(27) 


where  the  FIR  coefficients  matrix  is  A  =  [aoian...aii)T ,  and  an  auxiliary  quadratic  matrix  [( l  +  1)  x  (l  +  1)] 
is  given  by 


do 

di 

d2 

...  di 

di 

d'2 

d3 

...  di+1 

d2 

d3 

C?4 

...  di- 1_2 

.  di 

di- i-i 

di+ 2 

•  ••  d-2i 

N- 1 

where  the  generic  component  drn  =  im,  m  £  [0, 2Z] ,  is  determined  by  the  Bernoulli  polynomials  (Ap- 

i— o 

pendix  A).  An  analytic  solution  of  (27)  yields  the  generic  coefficients  for  (26) 


aji  =  (-1) 


jMU+ i)i 

|D| 


(29) 


in  which  |D|  and  are  the  determinant  and  minor  of  (26),  respectively.  Determined  aji  and  hi(n), 

the  unbiased  FIR  estimate  of  the  clock  ( K  +  1  —  Z)-th  state  is  given  to  be 


where 


N-l 

xK+i -i(n)  =  ^  hi(i)^K+i-i(n  -  i) 

i—0 


=  WfH  K+1-i(n), 


2 K+i-i{n i) 


€k+x -i(n) 
£,K+i-i{n  -  1) 

(. K+i-i(n  —  iV  +  1) 


(30) 

(31) 


(32) 


Estimate  Noise 

The  estimate  noise  variance  (18)  may  now  be  rewritten  as 


K 


var[A(n)]  =  ([r(n)r]T[r(n)r])  =  ^W^_iRfc+1(n)WJf.fe, 


k= 0 


R  i{n)  = 


Ri(n,n) 
Ri{n  -  1,  n) 


Ri(n,n~- 1) 
Ri{n  —  1,  n  —  1) 


...  Ri{n,  n  —  N  +  1) 

...  Ri(n  —  1,  n  —  N  +  1) 


(33) 


(34) 


Ri(n  —  N  +  l,n)  Ri(n  —  N  +  1,  n  —  1)  ...  Ri{n  —  N  +  1,  n  —  N  +  1) 

where  the  autocorrelation  matrix  R;(n)  is  specified  by  (33)  with  a  generic  component  Ri(i,j )  =  (vi(i)vi(j)), 
i,j&  [n,  n  —  N  +  1].  Accordingly,  the  estimate  variance  associated  with  the  estimate  (31)  calculates 

4'+i -iH  =  WfRK+1-t(n)W,.  (35) 


It  is  important  that,  by  large  N,  the  sawtooth  noise  becomes  delta-correlated.  This  degenerates  (34)  to  the 
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diagonal  form  with  the  components  Ri(i,i)  =  cr2;  and  (35)  is  written  as 

o-2K+1-l=a2v{K+^l)  WfWp  (36) 

where  ^(x+i-p  is  a  variance  of  the  noise  perturbing  the  (K  +  1  —  Z)-th  clock  state. 

Let  us  note  that,  in  different  algorithms,  the  components  of  the  noise  vector  v(n)  may  be  caused  by  different 
or  equal  physical  sources.  Thus,  they  may  demonstrate  different  powers  of  correlation  and  dependence. 


APPLICATIONS  TO  THE  CLOCK  TIE  POLYNOMIAL  MODEL 

In  applications,  K  in  (3)  is  identified  for  the  filter  memory  on  a  horizon  [0,7V  —  1])  by  the  clock  precision. 
Typically,  it  is  assumed  that  (3)  fits  cesium  and  hydrogen  maser  clocks  with  K  £  [0,1],  and  crystal  and 
rubidium  clocks  with  K  £  [1,2].  However,  K  =  3  may  be  required  for  low-precision  crystal  clocks.  Below 
we  give  the  unbiased  FIRs  for  all  these  cases. 


Low-Order  FIRs  for  the  Optimally  Unbiased  Filters 

Setting  l  =  0, 1,2,3  and  using  the  coefficient  dm  (Appendix),  we  first  calculate  the  FIR  coefficients  (29). 
Then  the  relevant  FIRs  (26)  may  readily  be  written  as 


ho(i )  =  jj, 


(37) 


hi(i )  = 


2(2 TV  -  1)  -  6i 


h2(i)  = 


7V(7V  + 1)  ’ 

3(3 TV2  -  37V  +  2)  -  18(27V  -  1  )i  +  30*2 


h3(i)  = 


7V(7V  +  1)(7V  +  2) 

8(27V3  -  3 TV2  +  7N  -  3)  -  20(67V2  -  67V  +  5 )i 


(38) 


(39) 


7V(7V  +  1)(7V  +  2)(7V- 
120(27V  -  l)i2  -  140«3 


3) 


'  7V(7V  +  1)(7 V  +  2) (TV  +  3) '  ^ 

The  constant  FIR  (37)  corresponds  to  simple  averaging  and  is  optimal  in  a  sense  of  minimum  produced 
noise.  This  FIR  is  practically  proved  to  be  reasonable  in  GPS-based  common  view  measurements  [5].  The 
linear  FIR  (38)  was  also  derived  in  [22]  by  using  linear  regression  to  compensate  a  bias  of  simple  averaging. 
This  filter  demonstrates  an  intermediate  error  between  the  2-D  and  3-D  Kalman  filters  in  application  to 


the  crystal  and  rubidium  clocks  [24].  Its  kernel  starts  with  a  maximum  h2(0)  = 


=  P  >« 


and  goes  to  a  minimum  h2(N  —  1)  =  —  ^(N+i) 


peculiarity  is  r  =  Jim  =  aU°ws  one  to  synthesize  a  FIR,  by  saving  r  =  —2  for  arbitrary 

TV.  It  may  be  shown  that  the  FIR  synthesized  in  such  a  way  is  equal  to  that  derived  in  [23]  for  the  1-step 
linear  prediction  on  a  horizon  [1,7V]. 


N»  1 


=  —  <  0,  having  zero  at  Hq  = 


N»1 

2N3~1  ■  Its  special 
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Noise  Power  Gains 

The  noise  power  gain  corresponding  to  the  Z-D  FIR  is  specified,  by  (36),  to  be  gi(N)  =  a = 
WfW i.  Its  values  associated  with  (37)-(40)  are  given  below,  respectively, 

9o(N)=±,  (41) 


92{N) 


9i(N) 


2{2N  —  1) 
N(N  +  1)  ’ 


3(3ZV2  -3 N+  2)  (IV2  +  3 AT  +  2) 
N{N  +  1)2(N  +  2)2 


8(2ZV3  -  3ZV2  +  7ZV  -  3) 
93^  ’  ~  N(N  +  1)(N  +  2)(N  +  3)  ’ 


(42) 

(43) 

(44) 


Relations  (41)-(44)  manifest  that  unbiasedness  is  achieved  at  an  increase  of  noise.  Indeed,  the  curves  for 
l  >  0  trace  above  the  lower  bound  1/y/N  associated  with  simple  averaging  (Z  =  0)  that  produces  minimum 
noise  (among  all  filters).  It  also  follows  that,  by  large  N  >  100,  the  noise  gain  is  performed  approximately 
by  \J 9i  (-ZV)  =  (Z  +  1)/VN  and,  thus,  traces  below  the  upper  bound;  that  is, 


Vgi(N)  < 


(1  +  1)/VN,N>{1  +  1)2 
1,N<  (Z  +  l)2. 


(45) 


ESTIMATING  THE  TIE  MODEL  WITH  A  SINGLE  GPS 
TIMING  RECEIVER 

We  now  consider  the  most  widely  used  practical  case  when  the  measurement  £i(n)  of  a  TIE  X\ (n)  is  obtained 
with  a  single  multichannel  GPS  timing  receiver  and  the  sawtooth  correction  is  not  used.  An  estimation 
algorithm  is  shown  in  Figure  2  for  the  K-D  polynomial  TIE  model  (3).  The  first  clock  state  estimate  Xi  (n) 
is  provided  with  hK{n)  by  (30)  at  a  horizon  [0,  NK  —  1],  The  observation  £2(n)  for  the  second  state  x2(n)  is 
then  formed,  using  (6),  as  a  discrete  time  derivative  of  xi (n).  Accordingly,  x2(n)  is  achieved  with  hx- i(u) 
at  a  horizon  [0,  Nk- i  —  1].  Inherently,  the  first  accurate  value  of  x2(n)  appears  at  (Nk  +  Nk-  \  —  2)th  point. 
Finally,  the  last  state  estimate  XK+i(n)  is  calculated  with  ho(n)  at  a  horizon  [0,ZVo  —  1]  using  ^+i  (n)  that 
is  formed  in  the  same  manner  as  t;2(n).  The  first  correct  value  of  Xk+i{ji)  appears  at  {Nk  +  Nk- i  +  ...  + 
N0  —  K  —  l)th  point. 

Below,  as  an  example  of  application,  we  use  this  algorithm  to  estimate  the  TIE  model  of  a  crystal  clock 
embedded  to  the  Stanford  Frequency  Counter  SR620.  The  measurement  is  done  with  SynPaQ  III  and  SR620 
for  A  =  1  s  (GPS-measurement).  Simultaneously,  to  get  a  reference  noiseless  trend,  the  TIE  of  the  same 
crystal  clock  is  measured,  by  SR625,  for  the  rubidium  clock  (Rb- measurement).  Initial  time  and  frequency 
shifts  between  two  measurements  are  then  eliminated  statistically  and  transition  to  A  =  10  s  is  provided 
by  the  data  thinning  in  time.  A  horizon  length  A)  for  each  FIR  filter  is  set  to  obtain  the  estimates  with 
minimum  MSEs. 
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Several-Hour  Measurement  and  Estimation 

In  this  experiment,  a  short-term  measurement  of  the  TIE  has  been  done  during  several  hours.  The  TIE  model 
was  then  identified  in  the  sense  of  a  minimum  MSE  to  be  quadratic.  Accordingly,  we  set  K  =  2,  estimate 
three  clock  states,  and  compare  the  results  for  the  Rb-measurement.  Figure  3  illustrates  the  estimates  X\(n) 
and  a:2(?i).  We  notice  that,  inherently,  a  several-hour  measurement  does  not  allow  for  a  proper  estimate  of 
the  third  clock  state  x3(n).  It  is  neatly  seen  (Figure  3a)  that  the  estimate  of  TIE  X\  (n)  follows  the  mean 
value  of  the  GPS-measurement  and  that  its  offset  from  the  Rb-measurement  is  caused  mostly  by  the  GPS 
time  uncertainty.  A  maximum  estimate  error  of  about  60  ns  has  appeared  in  the  time  span  between  the  8th 
and  9th  hours  when  a  time  shift  in  the  IPPS  signal  was  fixed.  Figure  3b  shows  that  the  estimate  a;2 (n)  of  the 
fractional  frequency  offset  fits  well  the  first  time  derivative  of  the  Rb-measurement.  However,  in  contrast 
to  (Figure  3a),  here  the  frequency  shift  of  about  3  x  10-11  has  occurred  in  the  span  between  the  7th  and 
8th  hours  and  no  appreciable  error  has  been  fixed  in  the  range  of  large  time  shift  (between  the  8th  and  9th 
hours) . 

Long-Term  Measurement  and  Estimation 

In  this  experiment,  we  watched  for  the  same  crystal  clock  over  about  2.5  days.  The  measurements  are 
inherently  fixed  oscillations  associated  with  a  day’s  variation  in  temperature  (Figure  4a).  We  notice  that, 
like  the  previous  case,  X\ (n)  fits  the  reference  Rb-measurement  with  the  error  caused  by  the  GPS  time 
uncertainty.  The  estimate  a;2  (n)  tracks  the  mean  value  of  the  first  time  derivative  of  Xi(n)  (Figure  46). 
Finally,  x2 (n)  fits  well  the  first  time  derivative  of  the  Rb-measurement  (Figure  4c).  And,  again,  a  2.5-day 
measurement  does  not  allow  for  estimating  x3{n)  with  a  sufficient  accuracy. 


CONCLUSION 

In  this  paper,  we  studied  an  optimally  unbiased  FIR  filter  of  the  TIE  polynomial  model  of  a  local  clock.  In 
contrast  to  the  standard  Kalman  filter,  the  proposed  solution  does  not  require  the  TIE  measured  process 
to  be  Gaussian  and  does  not  involve  any  knowledge  about  noise  in  the  algorithm.  Therefore,  the  algorithm 
seems  to  be  simpler  for  engineering  applications.  The  filter  is  asymptotically  optimal,  since  a  variance  of  the 
produced  noise  reduces  as  a  reciprocal  of  the  horizon  length  N.  Let  us  add  that  timekeeping  operates  typically 
with  large  horizons,  N  >>  1.  An  application  of  the  proposed  algorithm  to  the  GPS-based  measurement 
of  a  crystal  clock  showed  that  the  estimate  of  TIE  x\  fits  actual  values  with  an  uncertainty  of  the  GPS 
time.  Therefore,  such  a  filter  may  also  be  employed  as  an  estimator  of  the  GPS  time  uncertainty  in  hybrid 
structures.  The  estimate  of  the  fractional  frequency  offset  x2  fits  the  reference  Rb-measurement  with  high 
accuracy  that  is  limited  by  frequency  shifts  in  the  1  PPS  signal  of  the  GPS  receiver  and  in  the  reference 
signal.  Thus,  such  timing  receivers  may  efficiently  be  used  in  remote  measurements  of  frequency  offsets  of 
local  oscillators  instead  of  expensive  quantum  sources. 

APPENDIX:  THE  COEFFICIENTS  OF  MATRIX  (25) 

The  coefficients  for  (28)  are  calculated  by 

N~  1  j 

dm  =  V  [Bm+i(N)  -  Bm+1\ , 

TO  +  1 
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where  Bn(x)  is  the  Bernoulli  polynomial  and  Bn  =  Bn( 0)  is  the  Bernoulli  coefficient.  For  low  orders,  Bn(x ) 

may  be  found  in  reference  books.  For  high  orders,  the  following  recurrent  relation  is  valid: 

Bn(x)  =  n  j  Bn_i(x)dx  +  Bn. 

Several  low  order  coefficients  dm  are  given  below 

,_Ar  _  N(N  -  1)  J  N(N-1)(2N-1)  _  N2(N  —  l)2 

“o  — tv,  a i  — - - - ,  a2  — - - - ,  a3  — - - - , 

2  b  4 

N(N  —  1)(2N  —  1)(37V2  —  37V  —  1)  N2(N  -  1)2(2N2  -  21V  -  1) 

d4- - 30 - ’  d5  = - 12 - ’ 

N(N  -  1)  (21V  -  1)  (31V4  -  61V3  +  31V  +  1) 

— - . 

6  42 

For  a  large  horizon,  N  »  1,  the  coefficients  dm  may  be  calculated  by  dm\N>>1  =  . 
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V, 


Figure  1:  Typical  TIE  sawtooth  noise  structure  induced  by  a  GPS  timing  sensor  SynPaQ  III:  (a)  short-term, 
1  s  measurements  and  ( b )  long-term,  10  s  measurements  with  a  day’s  data  shifted  by  200  ns. 
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Figure  2:  Structure  of  the  unbiased  FIR  filtering  algorithm  for  the  A'-D  TIE  polynomial  model  observable 
with  a  single  GPS  timing  receiver. 
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Figure  3:  Short-time  measurement  and  estimation  of  the  crystal  clock  TIE  model  without  the  sawtooth 
correction:  (a)  GPS-measurement  of  the  TIE  (points),  Rb-measurement  of  the  TIE  (dotted),  and  quadratic 
unbiased  FIR  estimate  of  TIE  X\ (n)  (solid);  ( b )  Rb-measurement  of  the  fractional  frequency  offset  (points) 
and  its  GPS-based  linear  unbiased  FIR  estimate  x^  (n)  (solid)  (Note:  the  scale  here  in  10-11) 
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Figure  4:  Long-term  measurement  and  estimation  of  the  crystal  clock  TIE  without  the  sawtooth  correction: 
(a)  GPS-measurement  (points),  Rb-measurement  (dotted),  and  quadratic  unbiased  FIR  estimate  of  TIE 
X\ (n)  (solid);  ( b )  observation  £2{n)  (light)  and  linear  unbiased  FIR  estimate  of  the  fractional  frequency 
offset  x2 (n)  (bold)  (Note:  the  scale  here  in  10-10)  ;  (c)  Rb-measurement  (light)  and  estimate  x2(n)  (solid) 
(Note:  the  scale  here  in  10“ 10). 
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